OBSCURING FRACTION OF ACTIVE GALACTIC NUCLEI: IMPLICATIONS FROM RADIATION-DRIVEN FOUNTAIN MODELS

Active galactic nuclei (AGNs) such as quasars and Seyfert galaxies are important objects in the cosmic history. The enormous radiation energy of AGNs significantly influences the formation and evolution of galaxies and their environment. The evolution of AGNs in high redshift has been studied by wide-field optical and X-ray surveys, and these studies have recently become statistically more reliable based on large samples (e.g., Pâris et al. 2014; Ueda et al. 2014; Aird et al. 2015; Buchner et al. 2015). However, because AGNs are believed to be surrounded by an optically thick material for a fraction of the solid angles, we do not necessarily observe all the radiation from the nuclei. Because this obscuration may affect the apparent evolution of the luminosity function (Merloni & Heinz 2013), the fraction of absorbed objects is an essential quantity not only for interpreting observations (e.g., Lawrence & Elvis 1982; La Franca et al. 2005; Ballantyne et al. 2006; Hasinger 2008; Han et al. 2012) but also for theoretical models of AGNs evolution (e.g., Enoki et al. 2003, 2014; Fanidakis et al. 2011).

The fraction of obscured objects, hereafter called the obscured fraction f_obs, can be estimated by various observational methods. The nature of the obscuration should be different depending on the wavelength, central luminosity, and redshift, suggesting a complicated nature and evolution of the environment around AGNs (e.g., Lawrence & Elvis 2010; Merloni et al. 2014). In particular, f_obs inferred from optical and infrared observations is relatively high (f_obs ∼ 0.4–0.6; Lusso et al. 2013; Roseboom et al. 2013; Shao et al. 2013) and shows weak dependence on the luminosity, whereas X-ray samples suggest that f_obs decreases with increasing X-ray luminosity (Ueda et al. 2003, 2014; Hasinger 2008; Han et al. 2012; Merloni et al. 2014). This case may be understood using the receding torus model (Lawrence 1991; Simpson 2005), but recent surveys with large samples (XMM-Newton and Swift-BAT) suggested more complicated behavior in terms of L_X: f_obs is small for both low and high luminosities (i.e., there is a peak for a certain luminosity). The redshift evolution of f_obs is also suggested by different groups, although the results are not fully consistent each other (Brightman & Nandra 2011; Burlon et al. 2011; Ueda et al. 2014; Aird et al. 2015; Buchner et al. 2015).

The obscuring fraction as a function of the AGN luminosity and observed wavelength should be related to the morphology, internal structures, and size of the torus. If the obscured fraction of AGNs are related to the luminosity of AGNs, it would be natural to assume that such obscuring properties are directly caused by the radiative and mechanical feedback from the AGNs. As discussed by Elvis (2012), the weak correlation between the X-ray and optical obscuration might suggest that obscuration is occurred on different scales. The optical property of the X-ray-selected AGNs suggest that the radiation environment around the nucleus may determine the properties of the obscuration (Merloni et al. 2014).

The obscuration is not necessarily caused by static structures of the interstellar medium (ISM), such as the "donut-like" torus (Urry & Padovani 1995). Recently, Wada (2012) showed that radiation from the accretion disk drives vertical circulation of gas flows on the order of a few to tens of parsecs, forming thick torus-like structures. Based on this model, type-1 and -2 SEDs can be also consistently explained as a result of the outflows and backflows (Schartmann et al. 2014). This "radiation-driven fountain" has an advantage that it naturally explains a long-standing question about the torus paradigm: how does the torus maintain its thickness? At a much smaller scale, the outflows originating in the accretion disk could obscure the nucleus (Elvis 2000; Proga et al. 2000; Nomura et al. 2013). Moreover, Elitzur (2006) also claimed that the doughnut-like torus is not necessary, based on a similar picture.

In this paper, using three-dimensional radiation-hydrodynamic simulations, which are based on Wada (2012; hereafter Paper I) and Schartmann et al. (2014), we investigate how the obscured fraction depends on the properties of AGNs, i.e., luminosity and black hole (BH) mass (or Eddington ratio). In Section 2, the numerical methods and models are described, and the numerical results are shown in Section 3. The results are analytically discussed using a simple toy model in Section 4. The origin of the obscuring fraction suggested by recent surveys are discussed in Section 5, and the summary and remarks for the present study are given in Section 6. Here, we focus on the effect of the radiative feedback from the AGNs alone; possible synergetic effects with nuclear starbursts (e.g., Wada et al. 2002) will be discussed elsewhere.

2.1. Numerical Methods

2.2. Initial Conditions and Model Parameters

We follow the three-dimensional evolution of a rotating gas disk in a fixed spherical gravitational potential under the influence of radiation from the central source. To prepare quasi-steady initial conditions without radiative feedback, we first run models of axisymmetric and rotationally supported thin disks with uniform density profile. After the thin disks are dynamically settled, the radiation feedback is turned on.

Free parameters in the present simulations is the luminosity of the AGN (L_AGN) and the BH mass (M_BH). The Eddington ratios L_AGN/L_E, for the Eddington luminosity ${L}_{{\rm{E}}}\;\equiv 4\pi {{Gcm}}_{{\rm{p}}}{M}_{\mathrm{BH}}/{\sigma }_{T},$ are changed from 0.01 to 0.5 for three BH masses: ${M}_{\mathrm{BH},8}\equiv {M}_{\mathrm{BH}}/{10}^{8}{M}_{\odot }=0.13,1.3,$ and 13. The total AGN luminosity ranges 1.6 × 10⁴³–8.0 × 10⁴⁶ erg s⁻¹, which correspond to X-ray luminosities of ${L}_{{\rm{X}},44}\equiv {L}_{{\rm{X}}}/({10}^{44}\;\mathrm{erg}\;{{\rm{s}}}^{-1})$ = 0.01–8.35. We also prepare the initial gas disk with two different densities: ρ₀ = 1.2 and 3.6 × 10³ M_⊙ pc⁻³. Thus, the gas mass fraction to the BH mass ranges from 4 × 10⁻⁴ to 0.12. The AGN luminosity is assumed to be constant during calculations, which is typically a few million years.

We have modified some assumptions in Paper I to achieve more realistic models. As is in Paper I, the accretion disk itself is not solved in present simulations, and we assume an anisotropic radiation field for the central source. In contrast to the assumption in Paper I, i.e., ${F}_{\mathrm{AGN}}\propto \mathrm{cos}\theta$ for all wavelength, here the ultraviolet flux changes as ${F}_{\mathrm{UV}}(\theta )\propto \mathrm{cos}\theta (1+2\mathrm{cos}\theta )$ (Netzer 1987), where θ denotes the angle from the rotational axis (z-axis) in order to take into account the limb darkening effect of the thin accretion disk. As a result, the UV flux is about 20% larger for the z-axis and 25% smaller at θ ∼ 60° than the case in Paper I. Since the X-ray radiation would be originated in the hot coronal atmosphere of the accretion disk (Netzer 1987; Xu 2015), here we simply take spherical symmetry for the X-ray flux. In Paper I, L_X was fixed to 10% of the bolometric luminosity, but here we use luminosity-dependent correction for the X-rays (2–10 keV; Hopkins et al. 2007) and for the B-band (Marconi et al. 2004). The X-ray luminosity for M_BH = 1.3 × 10⁷ M_⊙ is about half of that assumed in Paper I. In Paper I, we also simply assumed that the hypothetical accretion disk is parallel to the galactic disk. However, in reality, any orientations of the accretion disk could be possible in dependent of the galactic disk (e.g., Kawaguchi & Mori 2010). Here we assumed that the direction of the emerging non-spherical radiation is inclined by 45° from the rotational axis of the circumnuclear gas disk in some models, and its three-dimensional gas dynamics is fully followed without assuming any symmetry.

The dust-to-gas ratio is assumed to be a widely used value, i.e., 0.01, instead of 1/30 in Paper I, although the actual ratio is not clear in the AGN environment.

As summarized in Table 1, we explored a total of 24 different sets of parameters (each model name represent parameters as described in the caption).

Table 1.  Model Parameters and Resultant Obscured Fraction

Models	${M}_{\mathrm{BH}}$ a	${\gamma }_{\mathrm{Edd}}$ b	${L}_{\mathrm{AGN},46}$ c	${L}_{{\rm{X}},44}$ d	${\phi }_{i}$ e	${\rho }_{0}^{{\rm{f}}}$	${f}_{\mathrm{obs},23}$ g	${f}_{\mathrm{obs},22}$ h
∗7L01	0.13	0.01	0.0016	0.01	0	0.40	0.11	0.33
7L05	0.13	0.05	0.008	0.04	0	0.40	0.14	0.37
7D30	0.13	0.3	0.016	0.18	0	1.22	0.17	0.28
7L20	0.13	0.2	0.032	0.13	0	0.40	0.13	0.24
7D20i	0.13	0.2	0.032	0.13	π/4	1.22	0.19	0.28
7L30	0.13	0.3	0.048	0.18	0	0.40	0.14	0.22
7D30	0.13	0.3	0.048	0.18	0	1.22	0.22	0.38
7D10	0.13	0.3	0.048	0.07	0	1.22	0.17	0.27
7L40	0.13	0.4	0.064	0.22	0	0.40	0.12	0.68
7D40i	0.13	0.4	0.064	0.22	π/4	1.22	0.16	0.63
7L40i	0.13	0.4	0.064	0.22	π/4	0.40	0.12	0.78
8L10	1.3	0.1	0.16	0.45	0	0.40	0.26	0.44
8L20	1.3	0.2	0.32	0.76	0	0.40	0.25	0.45
8D20	1.3	0.2	0.32	0.76	0	1.22	0.21	0.49
8L30	1.3	0.3	0.48	1.03	0	0.40	0.23	0.47
8D30	1.3	0.3	0.48	1.03	0	1.22	0.22	0.50
*8D30i	1.3	0.3	0.48	1.03	π/4	1.22	0.27	0.95
8L30i	1.3	0.3	0.48	1.03	π/4	0.40	0.28	0.79
8L40i	1.3	0.4	0.64	1.51	π/4	0.40	0.43	0.83
*8D50	1.3	0.5	0.8	1.51	0	1.22	0.21	0.66
8D50i	1.3	0.5	0.8	1.51	π/4	1.22	0.37	0.94
9D10	13	0.1	1.6	2.54	0	1.22	0.31	0.39
*9L40	13	0.4	6.4	7.08	0	0.40	0.25	0.36
9D50	13	0.5	8.0	8.35	0	1.22	0.36	0.48

Notes. Each model name represents "black hole mass" (7, 8, or 9) + density of the initial gas disk (L or D) + Eddington ratio (0.05–0.5) + inclination of the accretion disk (i: ϕ_i = π/4). Models with *in their names are shown in Figures 1 and 2.

^aBlack hole mass (10⁸M_⊙). ^bEddington ratio. ^cTotal AGN luminosity (10⁴⁶ erg s⁻¹). ^dX-ray luminosity (10⁴⁴ erg s⁻¹). ^eInclination angle of the accretion disk. ^fDensity of the initial disk (10³ M_⊙ pc⁻³). ^gObscuring fraction for column density N_H ≥ 10²³ cm⁻². ^hObscuring fraction for column density N_H ≥ 10²² cm⁻².

Download table as:  ASCIITypeset image

In this section, we derive the obscuring fraction for a given column density toward the nucleus using a quasi-steady three-dimensional density distribution under the influence of the radiation feedback. The results are summarized as a function of the X-ray luminosity.

Four representative models of the density structures of the gas at a quasi-stable state are shown in Figure 1. Figure 1(a) shows a typical radiation-driven fountain of the gas (Paper I), which consists of bipolar, non-steady, non-uniform outflows, failed winds (back flows to the disk), and a dense, thin disk. Physical origin of the fountain flows was partly discussed in Paper I, but here we add some analysis of the numerical results in Appendix A.1. See also the simple model in Section 4 for conditions to generate outflows. As discussed in Appendix A.2, the system reaches to a quasi-steady state of the gas circulation after about two rotational period at r = 10 pc, which is t ∼ 0.5 Myr for M_BH = 1.3 × 10³ M_⊙.

Zoom In Zoom Out Reset image size

Using the density distribution, the column density toward the nucleus as a function of the viewing angles (θ_v) and the azimuthal angles is plotted in Figure 2(a). In Figure 2, column densities for the models shown in Figure 1 are also plotted. Among these models, the column density tends to increase toward θ_v = 0 (edge-on) and reaches the maximum of N_H ≥ 10²⁵ cm⁻². This occurs because of the stable existence of the dense, thin gas disk for which the radiation from the central source does not affect the disk except for the central part. The bipolar outflows and the failed winds cover a fraction of solid angles for the central source in model 8D50 (Figure 1(a)). Therefore, the nucleus is obscured for $| {\theta }_{v}| \lesssim 90^\circ$ and $| {\theta }_{v}| \lesssim 20^\circ$ for N₂₂ ≡ 10²² cm⁻² and N₂₃ = 10 N₂₂, respectively. The density change for the azimuthal direction causes 0.5–1.0 dex difference in the column density for a given viewing angle, which causes a small ambiguity for the obscured fraction (f_obs) as a function of the column density.

Zoom In Zoom Out Reset image size

Figure 3 shows radial distribution of the cumulative column densities in a typical fountain model (8D50), from which we can know where the obscuration mostly occurs for a given line of sight. It shows that N_H ≳ 10²² cm⁻² at r = 16 pc (i.e., the edge of the computational box) is mostly attained the gas inside 5–10 pc, whose structures are resolved by 20–40 numerical grid points.

Zoom In Zoom Out Reset image size

In some low luminosity models (e.g., model 7L01: L_X,44 = 0.011 and M_BH,8 = 0.13) no steady outflows are observed (Figure 1(b)), and as a result, the opening angle of the obscuring material becomes large with f_obs(N₂₂) ≡ f_obs,22 = 0.33 (Figure 2(b)). On the other hand, the obscuring fraction also becomes small (f_obs,22 ≃ 0.4) for high luminosity with a massive BH (model 9L40: L_X,44 = 7.08 and M_BH,8 = 13), such as those shown in Figures 1(c) and 2(c). These high luminosity models do not also show steady bipolar outflows. This is because the acceleration due to the radiation pressure for dusty gas is not sufficiently large to reach the escape velocity for the luminous models, in which the dust sublimation radius moves outward ($\propto {L}_{\mathrm{AGN}}^{1/2}$). This is discussed using a kinematic model in Section 4.

As a result of the fully three-dimensional calculations, the central radiation source (the accretion disk) is not restricted to a parallel configuration with the rotational axis of the gas disk. In reality, the accretion disk should not necessarily be in the same plane of the circumnuclear gas disk. We found that if the axis of the central source is inclined by 45° from the rotational axis of the disk, and the Eddington ratio (γ_Edd) is relatively high (≳0.3), the majority of the solid angle for the AGN is covered by the fountain gas (Figures 1(d) and 2(d)).

The obscuring fraction of all 24 models are summarized in Table 1 and Figure 4, where f_obs for N₂₂ and N₂₃ are plotted for each model as a function of L_X,44. Figure 4 shows that nature of the obscuration depends on the luminosity of AGNs, and f_obs,22 is small in both low and high luminosity cases. Moreover, f_obs,22 shows a maximum value of 0.8–0.9 around L_X,44 ∼ 0.2–1.0. The obscuring fraction for N₂₃ gradually increases from f_obs,23 ≃ 0.1 to 0.3 in a range of L_X,44 = 0.01–10, which is almost the same as the initial conditions. This reflects that the scale height of the dense part of the disk is larger for more luminous models (i.e., more massive BHs). For a less dense column density, f_obs,22 increases from 0.3 to 0.5–0.6 toward L_X,44 ∼ 1.0, as a consequence of the formation of radiation-driven fountains. There is no significant difference on f_obs between models with different initial gas density ρ₀.

Zoom In Zoom Out Reset image size

We also find that f_obs,22 becomes quite large (∼0.7–0.9) for parameters that cause outflows (γ_edd ≳ 0.3 and M_BH,8 = 0.13 or 1.3) when the central accretion disk is inclined by 45° (Figures 1(d), 2(d), and open blue symbols in Figure 4). This occurs because the gas circulation between the outflows and the dense gas effectively covers a large fraction of the solid angle. On the other hand, for more luminous AGNs (L_X,44 > 2), the outflows are stalled; therefore, a relatively small obscured fraction is observed. For a given parameter set, note that f_obs should have an ambiguity of ∼10%, which comes from the azimuthal inhomogeneous structures of the fountains (see Figure 2).

Although the overall structures of N_H presented here is similar to those found in Paper I, there are some differences, because of the modified assumptions (Section 2.2). The radiative feedback assumed in this paper is less effective for generating outflows than those in Paper I for the same BH mass. Moreover the present models show smaller dispersion of N_H for a given viewing angle θ, especially for smaller θ. This is partly because the circumnuclear disks in the present model keep its geometrically thin shape in the computational box. On the other hand, the disks are followed by twice larger computational box in Paper I (i.e., 64³ pc). The disks tend to be gravitational unstable in the outer part (see discussion in Appendix A.2), where gravity of the BH is weaker, resulting in geometrically thick, inhomogeneous disk. This thick disks are not uniform for the azimuthal direction (see for example Figure 7 in Paper I), therefore the column density for a given viewing angle show dispersion. Although most obscuration is occurred in the central 10 pc, distributions and structures in the outer part of the disk on several tens parsecs scale may contribute to the observed column density, especially for the viewing angles close to the edge-on. Therefore, N_H could have larger dispersions, and as a result the obscuring fraction especially for N ≳ 10²³ cm⁻² presented here might be a "lower" limit with the contribution of the attenuation by the gas on larger scales.

As shown in Section 3 and Paper I, the three-dimensional gas dynamics under the effect of the non-spherical radiation are generally complicated. However, a simple kinematic model is helpful for understanding the behavior of the radiative feedback and the origin of the dependence of f_obs on L_AGN in Figure 4.

In the phenomena discussed here, the primary driving force producing the fountain of gas around the AGNs is the radiation pressure for the dusty gas. We assume that the accretion disk radiates most of its energy toward the direction of the rotational axis, not toward the plane of the gas disk (except for the case of the inclined accretion disk). On the other hand, the radiative heating more isotropically affects the surrounding gas through advection. The X-ray radiation from the AGN, which is assumed to be spherically symmetric in this study, heats the inner part of the thin gas disk, making it geometrically thick. Thus, we expect that the gas is radially accelerated by the radiation pressure with assistance from the X-ray heating. In the following discussion, the effect of the gaseous pressure in the outflowing gas and the self-gravity of the gas are ignored for simplicity. The outflowing gas is also dealt as a "parcel" or "shell," that is the attenuation of the central radiation by the outflows themselves are ignored.

Let us first assume that the X-ray heating causes the inner part of the gas disk to be a quasi-stable, uniform sphere with radius r₀, as schematically shown in Figure 5. In the gas sphere with r = r₀, the radiative cooling with cooling rate Λ_cool is assumed to be balanced with the X-ray heating rate: ${\rho }_{{\rm{g}}}^{2}{{\rm{\Lambda }}}_{\mathrm{cool}}\sim {L}_{{\rm{X}}}/(4\pi {r}_{0}^{3}/3),$ where ρ_g is the average gas density. Then, the radius r₀ within which the X-ray heating is effective is

$$\begin{eqnarray}&&{r}_{0}\sim {(3{L}_{{\rm{X}}}/4\pi {\rho }_{{\rm{g}}}^{2}{{\rm{\Lambda }}}_{\mathrm{cool}})}^{1/3}.\end{eqnarray} \tag{ 1 }$$

Zoom In Zoom Out Reset image size

Now, let us assume that the gas at r > r_d is accelerated due to the radiation pressure of the dust, where r_d is the dust sublimation radius.¹ Therefore, for the velocity v₀ at r = r₀, energy conservation requires

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}\left({v}_{0}^{2}-{v}_{\mathrm{max}}^{2}\right)={{GM}}_{\mathrm{BH}}\left(\displaystyle \frac{1}{{r}_{0}}-\displaystyle \frac{1}{{r}_{\mathrm{max}}}\right),\end{eqnarray} \tag{ 2 }$$

where r_max is the maximum radius that the launched gas reaches at maximum velocity v_max. For the dusty gas with opacity κ and dust-to-gas ratio γ_d, the gas can be accelerated toward v_max due to radiation pressure from r₀ to r_max:

$$\begin{eqnarray}&&{v}_{\mathrm{max}}^{2}\sim \displaystyle \frac{\kappa {\gamma }_{d}{L}_{\mathrm{AGN}}(\theta )}{4\pi c}\left(\displaystyle \frac{1}{{r}_{0}}-\displaystyle \frac{1}{{r}_{\mathrm{max}}}\right),\end{eqnarray} \tag{ 3 }$$

where L_AGN(θ) is the non-spherical radiation field of the AGN depending on the angle from the rotational axis θ. Similarly, v₀ is the velocity attained by to the radiation pressure from r_d to r₀, i.e.,

$$\begin{eqnarray}&&{v}_{0}^{2}\sim \displaystyle \frac{\kappa {\gamma }_{d}{L}_{\mathrm{AGN}}(\theta )}{4\pi c}\left(\displaystyle \frac{1}{{r}_{{\rm{d}}}}-\displaystyle \frac{1}{{r}_{0}}\right).\end{eqnarray} \tag{ 4 }$$

Using Equations (2)–(4), we obtain r_max as a function of L_AGN and M_BH:

$$\begin{eqnarray}&&{r}_{\mathrm{max}}({L}_{\mathrm{AGN}},{M}_{\mathrm{BH}})\;=\\ &&\quad \displaystyle \frac{{\hat{L}}_{\mathrm{AGN}}+{{GM}}_{\mathrm{BH}}}{{{GM}}_{\mathrm{BH}}/{r}_{0}-{\hat{L}}_{\mathrm{AGN}}(1/{r}_{{\rm{d}}}-2/{r}_{0})},\end{eqnarray} \tag{ 5 }$$

where ${\hat{L}}_{\mathrm{AGN}}\equiv \kappa {\gamma }_{{\rm{d}}}{L}_{\mathrm{AGN}}/(8\pi c).$ If ${r}_{\mathrm{max}}\to \infty ,$ the velocity of the dusty gas becomes larger than the escape velocity, and we expect that an outflow is formed.

For simplicity, we take ${L}_{\mathrm{AGN}}(\theta )=1/2{L}_{\mathrm{UV}}| \mathrm{cos}(\theta )| ,$ rather than ${L}_{\mathrm{AGN}}(\theta )\propto \mathrm{cos}\theta (1+2\mathrm{cos}\theta )$ as in the numerical simulations, where L_UV is the UV luminosity of the AGN. This difference is not essential for the following discussion. If ${{GM}}_{\mathrm{BH}}/{r}_{0}={\hat{L}}_{\mathrm{AGN}}(1/{r}_{{\rm{d}}}-2/{r}_{0})$ in Equation (5), then r_max tends toward infinity for the critical angle θ_crit, which is derived as

$$\begin{eqnarray}&&\mathrm{cos}{\theta }_{\mathrm{crit}}=\displaystyle \frac{{{GM}}_{\mathrm{BH}}}{{r}_{0}}\displaystyle \frac{16\pi c}{\kappa {\gamma }_{{\rm{d}}}{L}_{\mathrm{UV}}}{\left(\displaystyle \frac{1}{{r}_{{\rm{d}}}}-\displaystyle \frac{2}{{r}_{0}}\right)}^{-1}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&=\;\displaystyle \frac{{\sigma }_{{\rm{T}}}}{{m}_{{\rm{p}}}\kappa {\gamma }_{{\rm{d}}}\;{\gamma }_{\mathrm{Edd}}\;{r}_{0}}{\left(\displaystyle \frac{1}{{r}_{{\rm{d}}}}-\displaystyle \frac{2}{{r}_{0}}\right)}^{-1},\end{eqnarray} \tag{ 7 }$$

where γ_Edd is the Eddington ratio, σ_T is the Thomson cross-section, and m_p is the proton mass. We expect that, for θ < θ_crit, the radiation force is sufficiently large to allow the dusty gas to escape. Thus, we expect that, for this cone region, outflows of gas are formed. Otherwise, the gas launched from the inner region can only reach r = r_max(θ). We suppose that this gas eventually falls back toward the equatorial plane, which corresponds to the backflows seen in the hydrodynamic simulations (Figure 1(a) and Paper I).

Since both the equilibrium radius r₀ (Equation (1)) and the dust sublimation radius r_d depend on the luminosity of the AGN (e.g., ${r}_{{\rm{d}}}\propto {L}_{\mathrm{AGN}}^{0.5}$), r_max and θ_crit cannot simply be functions of the Eddington ratio γ_Edd, but they depend on both the AGN luminosity and the BH mass. For non-spherical radiation, r_max is also a function of the angle θ. In Figure 6, we plot r_max(θ) for five different parameter sets of L₄₆ ≡ 10⁴⁶ erg s⁻¹ and M_BH,8 ≡ 10⁸ M_⊙, showing the maximum radius of the launched gas. Here, we assume that ${r}_{{\rm{d}}}=1.3{L}_{46}^{1/2}$ pc (e.g., Lawrence 1991) and ${{\rm{\Lambda }}}_{\mathrm{cool}}={10}^{-22}/{n}_{{\rm{H}}}^{2}\;\mathrm{erg}\;{\mathrm{cm}}^{3}\;{{\rm{s}}}^{-2}.$ For (${L}_{46},{M}_{\mathrm{BH},8})=(0.02,\;0.1)$ and (0.2,1.0) (red and blue lines), r_max(θ) has finite values for some angles; thus, conical outflows are expected for θ < θ_crit. In the remaining three models with $({L}_{46},{M}_{\mathrm{BH},8})=(2.0,10.0),(10.0,10.0),$ and (0.01, 1.0), the accelerated gas cannot escape at any angle. In these cases, we suspect that steady outflows are not formed or they are stalled.

Zoom In Zoom Out Reset image size

In Figure 7, we plot θ_crit (Equation (7)) as a function of M_BH and L_X,44. Constant Eddington ratios are also shown with dotted lines. Here, we assume κ = 10³ cm² g⁻¹ and ρ_g = 300 M_⊙ pc⁻³. We also use the bolometric correction for X-rays (Hopkins et al. 2007) and for the B-band (Marconi et al. 2004). Figure 7 shows that there are regions where no outflows are formed, i.e., θ_crit has no value. For a given BH mass (M_BH ≲ 10⁹ M_⊙), there is a lower limit of the luminosity forming the outflows (i.e., π/2 > θ_crit > 0), under which the radiation pressure for the dusty gas is not sufficiently large to make the outflow velocity larger than the escape velocity. Interestingly, there is also an upper limit of the luminosity forming the outflows. This limit occurs because the dust sublimation radius r_d becomes close to the cooling radius r₀; as a result, the factor ${(1/{r}_{{\rm{d}}}-2/{r}_{0})}^{-1}$ is dominant in Equation (7).

Zoom In Zoom Out Reset image size

As we found in the radiation-driven fountain in the radiation-hydrodynamic simulations in Paper I and Section 2, part of the bipolar outflows driven by radiation from the AGNs falls back to the disk plane, forming a thick disk (the torus). The opening angle of the outflow can then be simply assumed to be close to the critical angle θ_crit derived in the above argument. Thus, for larger θ_crit, the obscuration fraction f_obs is expected to be smaller. On the other hand, if r_max is not infinity for all of the angles, the outflow is not formed. In this case, we expect that the central source is visible for most viewing angles.²

Based on these simple arguments, we can qualitatively understand the numerical results in Section 3. We found that the obscuration fraction for a given column density, especially for N₂₂, increases with the luminosity of the AGN, and turns to be a small value above L_X ≳ 10⁴⁵ erg s⁻¹. The increase of obscuration is consistent with the behavior of the opening angle as a function of the central luminosity (Figure 7). If the BH is too massive, e.g., M_BH > 10⁹M_⊙, however the outflow is not expected even for large luminosity. This situation occurs in the numerical simulations, and as a result, we observe a relatively small obscuration fraction for high luminosity AGNs.

In the simplest version of the receding torus model (Lawrence 1991; Hasinger 2008), a constant scale height is assumed for the gas disk. Then, the dust sublimation radius, which is a function of the central luminosity, is the only factor affecting the opening angle. However, our radiation hydrodynamics calculations and analytic argument provide a more dynamic picture. The radiation from the AGNs naturally drives outflows if conditions are satisfied, which affect the obscuring nature. The response of the dust to the central radiation is important, but this response would be more dynamically related with the obscuration than that assumed in the static receding torus model.

Note, however, that strong outflows are observed in some luminous quasars, which could originate in the line-driven wind from the accretion disk (Elvis 2000; Proga et al. 2000; Nomura et al. 2013). The line-driven outflows also contribute to the obscuration and could be related with ultra-fast outflows (UFOs).

5.1. Torus Structures Inferred from f_obs

5.2. Origin of Luminosity-dependent Obscuration in X-Rays

5.3. Implication to the Origin of the Redshift Evolution

The X-ray surveys discussed in the previous section also showed that obscured AGNs (Compton-thin AGNs) increase toward z = 3–5. For example, Buchner et al. (2015) showed that the fraction of the absorbed AGNs is ∼0.7 and ∼0.6 at L_X,44 ≃ 1 and 10, respectively, at z = 3 in contrast to ∼0.4–0.5 at z = 0.5. Aird et al. (2015) also suggested similar results and claimed that f_obs(L_X) differs depending on the luminosity at high-z. These studies also suggested that the fraction of Compton-thick objects is nearly constant (∼0.4 or 0.2) toward z = 3–5.

Figure 8 shows the opening angle of the bipolar outflows derived from the same analytical model as that shown in Figure 7. In this plot, θ_crit is plotted for three initial densities of the gas disk, i.e., ρ₃ = 1.0, 5.0, and 10, commonly showing that the opening angle becomes small as the central luminosity increases. Moreover, Figure 8 shows that, for a given luminosity, θ_crit is smaller for denser gas, i.e., more obscured. If the gas density in the central region is not too dense (otherwise, there is no solution for θ_crit, meaning that outflows are not formed), AGNs could be more obscured for more gas-rich environments. If high-z AGNs are associated with denser gas compared to local objects, this result is qualitatively consistent with the trend shown by the X-ray survey. Note again that this results is suggested only by the framework of the radiative feedback. In reality, we suspect that, in very gas-rich objects at high-z, nuclear starbursts should be initiated, which could be also play a role in the obscuration through supernova feedback (e.g., Levenson et al. 2001; Wada et al. 2002).

Zoom In Zoom Out Reset image size

Based on the radiation-driven fountain scenario (Wada 2012; Schartmann et al. 2014), which is the circulation of gas driven by the central radiation obscuring the nucleus, we investigated how the AGN properties, i.e., luminosity and BH mass, affect the obscuration due to the circumnuclear ISM, using fully three-dimensional radiation-hydrodynamic simulations. The obscuring fraction f_obs is derived from distribution of column density N_H toward the nucleus, using the numerical simulations.

We found that f_obs changes depending on N_H, the total luminosity of the AGN (L_AGN), and the BH mass (M_BH) (Section 3 and Figure 6). For a given L_AGN, f_obs is always smaller (∼0.2–0.3) for large column density (N_H > 10²³ cm⁻²) than that for a smaller column density N_H > 10²² cm⁻². This is a natural consequence of the fact that the dense gas surrounding the nuclei is mostly concentrated on a thin disk.

We found that f_obs for N_H = 10²² cm⁻² increases from ∼0.2 to ∼0.6 with the X-ray luminosity L_X in a range of L_X = 10^42–44 erg s⁻¹, and f_obs becomes small ∼0.4 at high luminosity (≥2 × 10⁴⁴ erg s⁻¹). At a moderate luminosity L_X ∼ 10⁴⁴ erg s⁻¹, a large fraction of the solid angles (f_obs > 0.5) can be obscured by the column density of N_H ≥ 10²² cm⁻², because of the outflows and thick disk caused by the back flows. The f_obs behavior has been simply demonstrated using a kinematic model (Section 4), in which the bipolar outflows are formed in specific ranges of L_AGN and M_BH (Figures 5 and 6), and their opening angles do not simply depend on the Eddington ratio.

The AGNs can be buried with f_obs > 0.7 for ${N}_{{\rm{H}}}\geqslant {10}^{22}$ cm⁻², when the central accretion disk is inclined by 45° from the rotational axis of the circumnuclear gas disk (Figures 1(d) and 2(d), and open symbols in Figure 3). Even in this configuration, the radiative feedback from the accretion disk cannot destroy a dense gas disk from a few to tens of parsecs. This causes non-spherical circulations of the gas (Figure 1(d)), and as a result, the majority of the solid angles are obscured by the gas with N_H ≥ 10²² cm⁻². The mass accretion can be possible through the dense gas disk, which is not completely destroyed, even in this case. This behavior could be a possible structure of "buried" AGNs.

The behavior of f_obs as a function of the X-ray luminosity is consistent with recent X-ray surveys (Burlon et al. 2011; Aird et al. 2015; Buchner et al. 2015) (Section 5.2). However, the luminosity for which f_obs reaches its maximum is a factor of ten smaller (i.e., ∼10⁴³ erg s⁻¹) than the observations for local AGNs, and the luminosity would more accurately fit those suggested in AGNs at z = 2–3 (Section 5.3).

Note that we did not consider the effect of stellar feedback around AGNs in the present simulations and analysis. Many AGNs are associated with circumnuclear starbursts (e.g., Esquej et al. 2014; Netzer 2015), and high-velocity outflows of molecular and ionized gas are often observed in quasars and ultra-luminous infrared galaxies. These activities could also change the structure of the circumnuclear disk, and as a result the obscuration of AGNs can be affected.

The author is grateful to R. Davies, M. Schartmann, Y. Ueda, and N. Kawakatsu for many helpful comments and discussions. The author also appreciates many valuable comments and suggestions by the anonymous referee. The numerical computations presented in this paper were performed on a Cray XC30 at the National Astronomical Observatory of Japan.

A.1. Mechanism of the Outflows

In this complementary section, we discuss dynamics and physical origin of the "radiation-driven fountain." One may notice in some models that vertical (z-direction) motions of the gas are driven (Figure 1). Since the radiation pressure on the dusty gas due to the central radiation source is solved only for the radial direction by the ray-tracing method in this paper, this is mainly due to the effect of the gas pressure gradient caused by the radiative and mechanical heating and also by the density gradient. In Figure 9, we show snapshots of density, temperature (T_g), Mach number (${\mathcal{M}}$) for vertical direction (i.e., ${\mathcal{M}}$ $\equiv \;{v}_{z}/{c}_{{\rm{s}}},$ where c_s is the sound speed), and pressure distributions on the x−z plane of model 8D50 (the inclination of the central accretion disk is zero). It shows that velocity of the gas in the bipolar outflows along the rotational axis (z-axis) is largest. This is natural because the flux due to the central source is largest for the direction. However, the directions of the outflows are not perfectly radial, especially near the central "core" (torus-like structure seen in the density map inside r ∼ 5 pc) and the edges of the conical outflows. As seen in the temperature map (Figure 9(b)), a hot core with T_g ∼ 10⁷ K is formed, which is directly heated by the X-ray from the central source. Its size roughly corresponds to the Compton radius ${r}_{{\rm{c}}}\equiv {{GM}}_{\mathrm{BH}}{m}_{{\rm{p}}}/{{kT}}_{{\rm{g}}}\;\simeq \;$ $4\;\mathrm{pc}$ ${({T}_{{\rm{g}}}/{10}^{7}\;{\rm{K}})}^{-1}$ $({M}_{\mathrm{BH}}/1.3\times {10}^{8}{M}_{\odot }).$ In general, thermally driven wind is expected where r > r_c (e.g., Begelman et al. 1983; Woods et al. 1996). This suggests that the gas motion just outside the hot core would be affected by the pressure gradient. However, as shown in Figures 9(c) and (d), the gases in the conical outflows are low temperature and highly supersonic (${\mathcal{M}}$ >10–20) at least for the vertical direction. This implies that the thermal pressure cannot keep driving the outflow motion, and the radiation pressure plays a key role as discussed in Section 4. This is an essential difference between the outflows observed here and the thermally driven "disk wind" from the accretion disk, which is a much smaller scale phenomenon and they are launched from the disk surface (Woods et al. 1996). The temperature of the gas disk in the present models are too low to generate the disk winds.

Zoom In Zoom Out Reset image size

It is also interesting that the "back flows" are generated outside the cone, whose direction is not necessarily radial. For example, the accreting gases change their direction near the disk plane, which is due to the thermal pressure as suggested by the small Mach number and also due to the relatively small radiation pressure near the disk plane. Self-gravity of the gas disk may also contribute to disturb the radial flows.

A.2. Quasi-stability of the System and Mass Accretion to the Central BH

Although the outflows driven by the radiation pressure is locally non-uniform and non-steady, we assume in Section 3 and Section 5 that the obscured nature is based on globally quasi-steady structures of the circumnuclear gas over ∼1 Myr. Figure 10 shows the gas mass in the outflows ($| z| \gt 5$ pc) and the central core region (r < 5 pc). It shows that the outflows are formed in the first ∼0.4 Myr, and a quasi-steady state is kept over 1 Myr. The gas mass in the central region is slightly decreasing with a rate of ∼0.2 M_⊙ yr⁻¹, which is mostly caused by the outflows (see also discussion below). The total gas mass in the system is about 6 × 10⁶ M_⊙, therefore the "life time" of the system would be ∼3 × 10⁷ years, if there is no mass supply from outside of the system.

Zoom In Zoom Out Reset image size

Figure 11 shows time evolution of the gas mass just outside the inner boundary (r ≤ 0.75 pc) of model 8D50 and 8D30i. It is significantly fluctuating (the time step is 50 years), and the distribution of the mass change rate ${{dM}}_{{\rm{g}}}/{dt}$ is roughly Gaussian as shown in the histogram. The root mean square (rms) of ${{dM}}_{{\rm{g}}}/{dt}$ is ∼0.3 M_⊙ yr⁻¹. Therefore if the short-period fluctuation of the gas mass caused by a temporal mass accretion toward the BH and immediate by supply from the disk, the rms accretion rate corresponds to about 10% Eddington accretion rate, which is comparable to the assumed luminosity of the central source. However, the average accretion rate over a longer period is much smaller; for example, for the "accretion phase" around t = 1 Myr (Figure 11), the accretion rate is about 2 × 10⁻³ M_⊙ yr⁻¹. Therefore, the present simulations suggest that the radiation-driven fountain, whose time scale is at least ∼Myr, cannot be sustained only by the mass accretion through the disk on a sub-pc scale.

Zoom In Zoom Out Reset image size

The mass accretion rate could be enhanced if the angular momentum transfer in the disk is more effective. This could be possible if the turbulent viscosity or gravitational torque due to the non-axisymmetric structures are enhanced due to gravitational instability (see for example Wada et al. 2002, and references therein). The Toomre's Q value, $Q\equiv \kappa {c}_{{\rm{s}}}/\pi G{{\rm{\Sigma }}}_{{\rm{g}}},$ where κ is epicyclic frequency and Σ_g is the surface density of the gas disk, is

$$\begin{eqnarray*}Q & \simeq & 0.5\left(\displaystyle \frac{{c}_{{\rm{s}}}}{1\;\mathrm{km}\;{{\rm{s}}}^{-1}}\right){\left(\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{g}}}}{1000\;{M}_{\odot }\;{\mathrm{pc}}^{-2}}\right)}^{-1}\\ & & \times {\left(\displaystyle \frac{{M}_{\mathrm{BH}}}{1.3\times {10}^{7}\;{M}_{\odot }}\right)}^{1/2}{\left(\displaystyle \frac{r}{10\;\mathrm{pc}}\right)}^{-3/2}.\end{eqnarray*}$$

Therefore the disk around M_BH = 1.3 × 10⁷ M_⊙ can be gravitationally unstable even for most regions, if the temperature of the gas is smaller than 100 K. The gas temperature inside r ≲ 5 pc is much larger than 10⁴ K due to the X-ray heating, which is assumed to be spherical, thus disk becomes stable and geometrically thick (see Figure 9(a) and (b)). On the other hand, outer disk (r > 10 pc) is expected to be gravitationally unstable especially for smaller BHs. The most unstable wavelength in the differentially rotating disk, ${\lambda }_{{\rm{J}}}\equiv 2{c}_{{\rm{s}}}^{2}/G{{\rm{\Sigma }}}_{{\rm{g}}}\simeq 1$ pc (T_g/100 K)(Σ_g/1000 M_⊙ )⁻¹, therefore the current spatial resolution, 0.25 pc, is merginal to resolve the gravitational instability of the coldest, thin gas disk (see also Figure 12). However, the "effective" sound velocity including the velocity dispersion of the self-gravitationally unstable disk can be much larger than 1 km s⁻¹ (see for example, Wada 2001; Wada et al. 2002), and also by perturbation due to the back flows. As a result, the outer disk cannot be very thin (see Figure 9) and the current spatial resolution is fine enough to resolve its vertical structures.

Zoom In Zoom Out Reset image size

In order to clarify how the spatial resolution affects the mass accretion toward the center through gravitational instability, we run two additional test calculations, as shown by Figure 12. Since the radiation pressure does almost no effect on the gas disk itself due to the angle dependence of the radiation field and to save computational time, here we follow initial evolution of the disk without the radiation feedback. We take a vertically smaller computational box, i.e., 32² × 8 pc, which is solved using 128² × 32 or 256² × 32 grid points. Although the dense disk is thinner in the model with 0.125 pc resolution, the spiral-like features originated in self-gravitational instability are resolved in both models. Evolution of the mass in the central region (r ≤ 0.75 pc) is also plotted for the two different resolution. We found that the average mass accretion rate is somewhat larger in the model with smaller resolution (0.125 pc), but the difference is not significant.

If the dense gas disk continues to a larger radius (e.g., r ∼ 100 pc) with similar surface density, the disk should be highly unstable, therefore larger mass inflow to the central region due to the non-axysymmetric structures or gravity-driven turbulence could be possible. This is an interesting issue to be explored using larger computational box and finer resolution in the near future, as well as the effect of supernova feedback (e.g., Wada & Norman 2002) with the radiative feedback.